Fatigue-life distribution (fatiguelife) — Birnbaum–Saunders model for fatigue failures

The fatigue-life distribution (also known as the Birnbaum–Saunders distribution) is a positive, right-skewed continuous distribution used to model time-to-failure under cyclic fatigue.

A key characterization is a normalizing transform to a standard normal:

\[ Z = \frac{1}{c}\left(\sqrt{\frac{X}{\beta}} - \sqrt{\frac{\beta}{X}}\right) \sim \mathcal N(0,1), \]

equivalently, a constructive representation:

\[ X = \beta\left(\frac{cZ}{2} + \sqrt{\left(\frac{cZ}{2}\right)^2 + 1}\right)^2,\qquad Z\sim\mathcal N(0,1). \]

---

Learning goals

• understand what fatiguelife models and when it is used

• write down the PDF/CDF and connect them to a standard normal

• compute moments (mean/variance/skewness/kurtosis) and interpret parameters

• derive expectation, variance, and the (log-)likelihood

• sample from the distribution using NumPy only

• visualize PDF/CDF and Monte Carlo simulations

• use scipy.stats.fatiguelife for pdf, cdf, rvs, fit

import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
from scipy import stats, special
from scipy.stats import fatiguelife as fatiguelife_sp

# Plotly rendering (CKC convention)
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

# Reproducibility
rng = np.random.default_rng(7)
np.set_printoptions(precision=6, suppress=True)

SQRT2 = np.sqrt(2.0)
SQRT2PI = np.sqrt(2.0 * np.pi)

1) Title & classification

• Name: fatiguelife

• Type: continuous distribution

• Standard support: \(x \in (0, \infty)\)

• SciPy parameter space (scipy.stats.fatiguelife(c, loc=0, scale=1)):

  • shape: \(c > 0\) (often written \(\alpha\))

  • location: \(\mathrm{loc} \in \mathbb{R}\)

  • scale: \(\mathrm{scale} > 0\) (often written \(\beta\))

• Support with loc/scale: \(x \in (\mathrm{loc}, \infty)\)

We’ll use \((c,\beta)\) for the Birnbaum–Saunders form and note that in SciPy: \(\beta = \mathrm{scale}\).

2) Intuition & motivation

What it models

The fatigue-life distribution was introduced to model lifetimes of materials under repeated stress cycles. Intuitively, a material accumulates microscopic damage over many cycles; failure happens when accumulated damage crosses a threshold. This leads to positive support and typically right-skewed lifetimes.

Typical real-world use cases

• Reliability engineering: time-to-failure (or number of cycles to failure) under fatigue

• Materials science: crack initiation/propagation times

• Survival analysis: positive, skewed durations where lognormal is a competitor

Relations to other distributions

• Normal connection: the transform $\(Z = \frac{1}{c}\left(\sqrt{\frac{X}{\beta}} - \sqrt{\frac{\beta}{X}}\right)\)$ is standard normal.

• Lognormal-like behavior: for small \(c\), the distribution concentrates near \(\beta\) and can look approximately lognormal.

• Reciprocal symmetry: if \(X\sim\mathrm{BS}(c,\beta)\) then \(\beta^2/X\) has the same distribution. A nice consequence is that the median is exactly \(\beta\).

3) Formal definition

We write \(X \sim \mathrm{BS}(c,\beta)\) (Birnbaum–Saunders) with \(c>0\) and \(\beta>0\). Let \(\Phi\) and \(\varphi\) be the standard normal CDF and PDF.

CDF

For \(x>0\):

\[ F(x\mid c,\beta) = \Phi\!\left(\frac{1}{c}\left(\sqrt{\frac{x}{\beta}} - \sqrt{\frac{\beta}{x}}\right)\right). \]

PDF

For \(x>0\):

\[ f(x\mid c,\beta) = \varphi\!\left(\frac{1}{c}\left(\sqrt{\frac{x}{\beta}} - \sqrt{\frac{\beta}{x}}\right)\right) \cdot \frac{1}{2cx}\left(\sqrt{\frac{x}{\beta}} + \sqrt{\frac{\beta}{x}}\right). \]

Location–scale form (SciPy)

SciPy uses the standard rv_continuous transformation:

\[X = \mathrm{loc} + \mathrm{scale}\cdot Y,\quad Y\sim \mathrm{BS}(c,1).\]

So for \(x>\mathrm{loc}\), with \(y=(x-\mathrm{loc})/\mathrm{scale}\):

\[ f_X(x) = \frac{1}{\mathrm{scale}}\,f_Y(y),\qquad F_X(x)=F_Y(y). \]

def _validate_params(c: float, scale: float) -> None:
    if not (np.isfinite(c) and c > 0):
        raise ValueError("c must be finite and > 0")
    if not (np.isfinite(scale) and scale > 0):
        raise ValueError("scale must be finite and > 0")


def fatiguelife_pdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    '''Fatigue-life PDF in SciPy's loc/scale parameterization.

    Notes
    -----
    Uses the standard form with β=1 on y=(x-loc)/scale and applies the Jacobian 1/scale.
    '''
    _validate_params(float(c), float(scale))

    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    out = np.zeros_like(y, dtype=float)
    mask = y > 0
    if not np.any(mask):
        return out

    y_m = y[mask]
    s = np.sqrt(y_m)
    z = (s - 1.0 / s) / c

    phi = np.exp(-0.5 * z**2) / SQRT2PI
    jac = (s + 1.0 / s) / (2.0 * c * y_m)
    out[mask] = (phi * jac) / scale
    return out


def fatiguelife_logpdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    '''Log-PDF (useful for numerical stability).'''
    _validate_params(float(c), float(scale))

    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    out = np.full_like(y, fill_value=-np.inf, dtype=float)
    mask = y > 0
    if not np.any(mask):
        return out

    y_m = y[mask]
    s = np.sqrt(y_m)
    z = (s - 1.0 / s) / c

    log_phi = -0.5 * z**2 - np.log(SQRT2PI)
    log_jac = np.log(s + 1.0 / s) - np.log(2.0 * c * y_m)
    out[mask] = log_phi + log_jac - np.log(scale)
    return out


def fatiguelife_cdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    '''Fatigue-life CDF in SciPy's loc/scale parameterization.'''
    _validate_params(float(c), float(scale))

    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    out = np.zeros_like(y, dtype=float)
    mask = y > 0
    if not np.any(mask):
        return out

    y_m = y[mask]
    s = np.sqrt(y_m)
    z = (s - 1.0 / s) / c

    out[mask] = special.ndtr(z)
    return out


def fatiguelife_sf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    '''Survival function 1-CDF, computed stably via log_ndtr when possible.'''
    _validate_params(float(c), float(scale))

    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    out = np.ones_like(y, dtype=float)
    mask = y > 0
    if not np.any(mask):
        return out

    y_m = y[mask]
    s = np.sqrt(y_m)
    z = (s - 1.0 / s) / c

    # For large z, 1 - ndtr(z) loses precision; use log_ndtr(-z).
    out[mask] = np.exp(special.log_ndtr(-z))
    return out

4) Moments & properties

For \(X\sim\mathrm{BS}(c,\beta)\) with \(c>0\) and \(\beta>0\):

Mean, variance, skewness, kurtosis

• Mean $\(\mathbb E[X] = \beta\left(1 + \frac{c^2}{2}\right).\)$

• Variance $\(\mathrm{Var}(X)=\beta^2\,c^2\left(1+\frac{5c^2}{4}\right).\)$

• Skewness $\(\gamma_1 = \frac{4c\,(11c^2+6)}{(5c^2+4)^{3/2}}.\)$

• Excess kurtosis $\(\gamma_2 = \frac{6c^2\,(93c^2+40)}{(5c^2+4)^2}.\)$

(Scale and location do not change skewness/kurtosis.)

MGF and characteristic function

A simple closed form is not available. A useful representation uses the normal construction. Let \(Z\sim\mathcal N(0,1)\) and define

\[Y = \left(\frac{cZ}{2} + \sqrt{\left(\frac{cZ}{2}\right)^2 + 1}\right)^2.\]

Then for SciPy’s parameters \(X=\mathrm{loc}+\mathrm{scale}\,Y\):

• MGF (when it exists): $\(M_X(t)=\mathbb E[e^{tX}] = e^{t\,\mathrm{loc}}\,\mathbb E\left[e^{t\,\mathrm{scale}\,Y}\right].\)\( Because \)Y\( grows like \)\approx c^2Z^2\( for large positive \)Z\(, the MGF exists at least for \)\(t < \frac{1}{2\,\mathrm{scale}\,c^2}.\)$

• Characteristic function (always exists): $\(\varphi_X(\omega)=\mathbb E[e^{i\omega X}].\)$

Entropy

The differential entropy \(H(X)=-\mathbb E[\log f(X)]\) has no common elementary closed form. In practice you compute it numerically (e.g. via scipy.stats.fatiguelife.entropy or Monte Carlo).

def fatiguelife_moments(c: float, loc: float = 0.0, scale: float = 1.0) -> dict:
    '''Closed-form mean/var/skew/excess-kurtosis for SciPy's parameters.'''
    _validate_params(float(c), float(scale))

    mean = loc + scale * (1.0 + 0.5 * c**2)
    var = (scale**2) * (c**2) * (1.0 + 1.25 * c**2)
    skew = 4.0 * c * (11.0 * c**2 + 6.0) / ((5.0 * c**2 + 4.0) ** 1.5)
    excess_kurt = 6.0 * c**2 * (93.0 * c**2 + 40.0) / ((5.0 * c**2 + 4.0) ** 2)
    median = loc + scale * 1.0  # median of BS(c,1) is 1

    return {
        'mean': mean,
        'var': var,
        'skew': skew,
        'excess_kurt': excess_kurt,
        'median': median,
    }


def normal_expectation_gauss_hermite(func, n: int = 60) -> float:
    '''Approximate E[func(Z)] for Z~N(0,1) via Gauss–Hermite quadrature (NumPy-only).'''
    xs, ws = np.polynomial.hermite.hermgauss(n)
    # For standard normal: E[g(Z)] = 1/sqrt(pi) * sum w_i * g(sqrt(2) * x_i)
    zs = SQRT2 * xs
    return float(np.sum(ws * func(zs)) / np.sqrt(np.pi))


def fatiguelife_from_normal(z, c: float, loc: float = 0.0, scale: float = 1.0):
    '''Deterministic transform mapping N(0,1) samples -> fatigue-life samples.'''
    _validate_params(float(c), float(scale))

    z = np.asarray(z, dtype=float)
    delta = 0.5 * c * z
    y = delta + np.sqrt(delta**2 + 1.0)
    return loc + scale * (y**2)


def fatiguelife_mgf_gh(t: float, c: float, loc: float = 0.0, scale: float = 1.0, n: int = 80) -> float:
    '''MGF via Gauss–Hermite quadrature (best for moderate t).'''

    def g(z):
        x = fatiguelife_from_normal(z, c=c, loc=0.0, scale=scale)
        return np.exp(t * x)

    return np.exp(t * loc) * normal_expectation_gauss_hermite(g, n=n)


def fatiguelife_cf_gh(omega: float, c: float, loc: float = 0.0, scale: float = 1.0, n: int = 80) -> complex:
    '''Characteristic function via Gauss–Hermite quadrature.'''

    def g(z):
        x = fatiguelife_from_normal(z, c=c, loc=0.0, scale=scale)
        return np.exp(1j * omega * x)

    return np.exp(1j * omega * loc) * normal_expectation_gauss_hermite(g, n=n)


params = dict(c=1.0, loc=0.0, scale=2.0)
closed = fatiguelife_moments(**params)
scipy_stats = fatiguelife_sp.stats(params['c'], loc=params['loc'], scale=params['scale'], moments='mvsk')
closed, scipy_stats

({'mean': 3.0,
  'var': 9.0,
  'skew': 2.5185185185185186,
  'excess_kurt': 9.851851851851851,
  'median': 2.0},
 (3.0, 9.0, 2.5185185185185186, 9.851851851851851))

5) Parameter interpretation

Shape parameter \(c\) (= c)

• Controls spread, right tail weight, and skewness.

• As \(c\to 0\), the distribution collapses toward its median (\(\beta\) in the Birnbaum–Saunders form).

• Larger \(c\) produces heavier right tails and larger skewness/kurtosis.

Scale parameter \(\beta\) (= scale)

• Multiplies the variable: if \(Y\sim\mathrm{BS}(c,1)\) then \(X=\beta Y\sim\mathrm{BS}(c,\beta)\).

• The median equals \(\beta\).

Location parameter loc

• Shifts the support to \((\mathrm{loc},\infty)\).

• In reliability contexts, you sometimes interpret loc as a minimum/threshold lifetime.

x = np.linspace(1e-3, 12, 800)

param_sets = [
    (0.2, 0.0, 1.0, 'c=0.2, scale=1'),
    (0.5, 0.0, 1.0, 'c=0.5, scale=1'),
    (1.0, 0.0, 1.0, 'c=1.0, scale=1'),
    (2.0, 0.0, 1.0, 'c=2.0, scale=1'),
]

fig = go.Figure()
for c, loc, scale, label in param_sets:
    fig.add_trace(go.Scatter(x=x, y=fatiguelife_pdf(x, c=c, loc=loc, scale=scale), mode='lines', name=label))

fig.update_layout(
    title='Fatigue-life PDF: effect of the shape parameter c',
    xaxis_title='x',
    yaxis_title='density',
    width=900,
    height=420,
)
fig

6) Derivations

Expectation

Start from the constructive representation with \(Z\sim\mathcal N(0,1)\) and \(\delta = \tfrac{cZ}{2}\):

\[ Y = \left(\delta + \sqrt{\delta^2+1}\right)^2 = 1 + 2\delta^2 + 2\delta\sqrt{\delta^2+1}. \]

Take expectations. The last term is an odd function of \(Z\) (because \(\delta\) is odd and \(\sqrt{\delta^2+1}\) is even), so it integrates to zero under a symmetric normal.

\[ \mathbb E[Y] = 1 + 2\mathbb E[\delta^2] = 1 + 2\cdot\frac{c^2}{4}\mathbb E[Z^2] = 1 + \frac{c^2}{2}. \]

Finally, for SciPy parameters \(X=\mathrm{loc}+\mathrm{scale}\,Y\):

\[\mathbb E[X] = \mathrm{loc} + \mathrm{scale}\left(1 + \frac{c^2}{2}\right).\]

Variance

Write \(Y=A+B\) with

• \(A = 1+2\delta^2\) (even in \(Z\))

• \(B = 2\delta\sqrt{\delta^2+1}\) (odd in \(Z\))

Then

\[Y^2 = A^2 + 2AB + B^2,\]

and \(\mathbb E[AB]=0\) because \(AB\) is odd in \(Z\). Using \(\mathbb E[Z^2]=1\) and \(\mathbb E[Z^4]=3\), you get

\[ \mathbb E[Y^2] = 1 + 2c^2 + \frac{3c^4}{2}. \]

So

\[ \mathrm{Var}(Y)=\mathbb E[Y^2] - (\mathbb E[Y])^2 = c^2\left(1+\frac{5c^2}{4}\right). \]

Scaling gives

\[\mathrm{Var}(X)=\mathrm{scale}^2\,\mathrm{Var}(Y).\]

Likelihood (i.i.d. sample)

Given data \(x_1,\dots,x_n\) and parameters \((c,\mathrm{loc},\mathrm{scale})\), define \(y_i=(x_i-\mathrm{loc})/\mathrm{scale}\). Validity requires \(c>0\), \(\mathrm{scale}>0\), and all \(y_i>0\).

For each \(i\):

\[ \log f(x_i) = -\log(\mathrm{scale}) -\log(2c) - \log(y_i) + \log\left(\sqrt{y_i}+\frac{1}{\sqrt{y_i}}\right) -\tfrac12 z_i^2 - \tfrac12\log(2\pi), \]

where \(z_i=\tfrac{1}{c}(\sqrt{y_i} - 1/\sqrt{y_i})\). The total log-likelihood is the sum over \(i\).

def fatiguelife_nll(params, data):
    '''Negative log-likelihood for (c, loc, scale).'''
    c, loc, scale = params
    if not (np.isfinite(c) and c > 0 and np.isfinite(scale) and scale > 0 and np.isfinite(loc)):
        return np.inf

    y = (data - loc) / scale
    if np.any(y <= 0):
        return np.inf

    return -float(np.sum(fatiguelife_logpdf(data, c=c, loc=loc, scale=scale)))


# Quick check: compare our logpdf to SciPy's for random x>0
c0, loc0, scale0 = 0.8, 0.0, 1.7
x_test = rng.uniform(1e-4, 8, size=10_000)
err = np.max(
    np.abs(
        fatiguelife_logpdf(x_test, c=c0, loc=loc0, scale=scale0)
        - fatiguelife_sp.logpdf(x_test, c0, loc=loc0, scale=scale0)
    )
)
err

1.1368683772161603e-12

7) Sampling & simulation (NumPy only)

The normal-construction gives a very clean sampler.

Algorithm

To sample \(X\sim\mathrm{BS}(c,\beta)\) (with optional loc):

1. Sample \(Z\sim\mathcal N(0,1)\).

2. Compute \(\delta = \tfrac{cZ}{2}\).

3. Compute \(Y = \delta + \sqrt{\delta^2+1}\).

4. Return \(X = \mathrm{loc} + \beta\,Y^2\).

This is vectorized and uses only numpy.random + elementary operations.

def fatiguelife_rvs_numpy(c: float, loc: float = 0.0, scale: float = 1.0, size=1, rng=None):
    '''Sample from the fatigue-life distribution using NumPy only.'''
    _validate_params(float(c), float(scale))

    if rng is None:
        rng = np.random.default_rng()

    z = rng.standard_normal(size=size)
    return fatiguelife_from_normal(z, c=c, loc=loc, scale=scale)


samples_np = fatiguelife_rvs_numpy(1.0, loc=0.0, scale=1.5, size=5, rng=rng)
samples_np

array([2.243943, 1.852748, 0.772983, 7.725928, 0.481773])

8) Visualization

We’ll compare:

• the theoretical PDF/CDF

• Monte Carlo samples (histogram + empirical CDF)

c0, loc0, scale0 = 1.0, 0.0, 2.0
n = 80_000
samps = fatiguelife_rvs_numpy(c0, loc=loc0, scale=scale0, size=n, rng=rng)

x_grid = np.linspace(1e-3, np.quantile(samps, 0.995), 600)

# Empirical CDF
xs_sorted = np.sort(samps)
ecdf = np.arange(1, n + 1) / n

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=('PDF + Monte Carlo histogram', 'CDF + empirical CDF'),
)

# PDF panel
hist_y, hist_x = np.histogram(samps, bins=120, density=True)
hist_xc = 0.5 * (hist_x[1:] + hist_x[:-1])
fig.add_trace(go.Bar(x=hist_xc, y=hist_y, name='MC histogram', opacity=0.45), row=1, col=1)
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fatiguelife_pdf(x_grid, c=c0, loc=loc0, scale=scale0),
        mode='lines',
        name='theoretical pdf',
        line=dict(width=3),
    ),
    row=1,
    col=1,
)

# CDF panel
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fatiguelife_cdf(x_grid, c=c0, loc=loc0, scale=scale0),
        mode='lines',
        name='theoretical cdf',
        line=dict(width=3),
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=xs_sorted[::40],
        y=ecdf[::40],
        mode='markers',
        name='empirical cdf',
        marker=dict(size=4),
    ),
    row=1,
    col=2,
)

fig.update_layout(width=1000, height=420, legend=dict(orientation='h', y=1.15))
fig.update_xaxes(title_text='x', row=1, col=1)
fig.update_yaxes(title_text='density', row=1, col=1)
fig.update_xaxes(title_text='x', row=1, col=2)
fig.update_yaxes(title_text='F(x)', row=1, col=2)
fig

9) SciPy integration (scipy.stats.fatiguelife)

SciPy’s object supports the standard continuous-distribution API:

• pdf, logpdf, cdf, sf

• rvs for sampling

• fit for parameter estimation

• stats, moment, entropy, …

The call signature is:

scipy.stats.fatiguelife(c, loc=0, scale=1)

where c is the shape parameter and scale corresponds to \(\beta\).

# Sampling via SciPy
samps_sp = fatiguelife_sp.rvs(c0, loc=loc0, scale=scale0, size=n, random_state=rng)

# Fit parameters (returns: c_hat, loc_hat, scale_hat)
fit_c, fit_loc, fit_scale = fatiguelife_sp.fit(samps_sp)
fit_c, fit_loc, fit_scale

(0.9997818369711162, 0.001552102122452088, 2.0038068173424053)

# Compare fitted PDF to the true PDF
x_grid = np.linspace(1e-3, np.quantile(samps_sp, 0.995), 600)

fig = go.Figure()
fig.add_trace(go.Histogram(x=samps_sp, nbinsx=120, histnorm='probability density', name='data', opacity=0.45))
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fatiguelife_sp.pdf(x_grid, c0, loc=loc0, scale=scale0),
        mode='lines',
        name='true pdf',
        line=dict(width=3),
    )
)
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fatiguelife_sp.pdf(x_grid, fit_c, loc=fit_loc, scale=fit_scale),
        mode='lines',
        name='fitted pdf',
        line=dict(width=3, dash='dash'),
    )
)
fig.update_layout(
    title='SciPy fit: true vs fitted PDF',
    xaxis_title='x',
    yaxis_title='density',
    width=900,
    height=420,
)
fig

10) Statistical use cases

Hypothesis testing (goodness-of-fit)

A common approach is the probability integral transform (PIT). If the model is correct and parameters are known, then

\[U_i = F(X_i) \sim \mathrm{Uniform}(0,1).\]

So we can test uniformity of \(U_i\) using a KS test.

Important caveat: if you fit parameters on the same data and then run KS, the p-values are not exact. A more principled approach is a parametric bootstrap.

Bayesian modeling

For positive parameters (\(c>0\), \(\beta>0\)) it’s common to model them on a log-scale, e.g.

\[\log c \sim \mathcal N(\mu_c,\sigma_c^2),\qquad \log \beta \sim \mathcal N(\mu_\beta,\sigma_\beta^2).\]

Then the posterior is proportional to

\[p(c,\beta\mid x) \propto p(x\mid c,\beta)\,p(c)\,p(\beta).\]

Below we’ll do a simple grid posterior demonstration (no external Bayesian libraries).

Generative modeling

fatiguelife is a convenient generator for positive, skewed data. For example, it can model synthetic failure times in simulations, or serve as a base distribution in mixture/latent-variable models.

# Goodness-of-fit via PIT + KS (illustrative)
# We'll use the fitted SciPy params from above.
U = fatiguelife_sp.cdf(samps_sp, fit_c, loc=fit_loc, scale=fit_scale)
ks = stats.kstest(U, 'uniform')
ks

KstestResult(statistic=0.0016533666280400539, pvalue=0.9807346694212495, statistic_location=0.7203466333719599, statistic_sign=1)

# Simple Bayesian grid posterior for (c, scale) with loc fixed at 0
# (This is a pedagogical demo; for serious work use MCMC/VI.)

loc_fixed = 0.0

# Keep a smaller dataset for speed
data = samps_sp[:2000]

# Grid over log-parameters
logc_grid = np.linspace(np.log(0.2), np.log(3.0), 120)
logs_grid = np.linspace(np.log(0.3), np.log(6.0), 120)

# Log-likelihood on the grid (vectorized over scale for each c)
ll = np.zeros((logc_grid.size, logs_grid.size), dtype=float)
scale_grid = np.exp(logs_grid)
for i, c in enumerate(np.exp(logc_grid)):
    ll[i, :] = np.sum(fatiguelife_sp.logpdf(data[:, None], c, loc=loc_fixed, scale=scale_grid[None, :]), axis=0)

# Priors: log c ~ N(0, 1^2), log scale ~ N(0, 1^2)
log_prior = -0.5 * (logc_grid[:, None] ** 2 + logs_grid[None, :] ** 2)

log_post_unnorm = ll + log_prior
log_post = log_post_unnorm - np.max(log_post_unnorm)
post = np.exp(log_post)
post /= np.sum(post)

# Posterior summaries
c_mean = float(np.sum(np.exp(logc_grid)[:, None] * post))
scale_mean = float(np.sum(np.exp(logs_grid)[None, :] * post))
(c_mean, scale_mean)

(1.0045700659359345, 1.9890536073446352)

# Visualize the posterior over (c, scale)
fig = go.Figure(
    data=go.Heatmap(
        x=np.exp(logs_grid),
        y=np.exp(logc_grid),
        z=post,
        colorscale='Viridis',
        colorbar=dict(title='posterior mass'),
    )
)

fig.update_layout(
    title='Grid posterior p(c, scale | data) with log-normal priors (loc fixed at 0)',
    xaxis_title='scale (β)',
    yaxis_title='c',
    width=850,
    height=520,
)

# Add markers for truth and posterior mean
fig.add_trace(go.Scatter(x=[scale0], y=[c0], mode='markers', name='true', marker=dict(size=10, symbol='x')))
fig.add_trace(
    go.Scatter(
        x=[scale_mean],
        y=[c_mean],
        mode='markers',
        name='posterior mean',
        marker=dict(size=10, symbol='circle-open'),
    )
)

fig

11) Pitfalls

• Invalid parameters: must have c > 0, scale > 0. With loc, all observations must satisfy \(x>\mathrm{loc}\).

• Parameterization confusion: many texts use \((\alpha,\beta)\) for (shape, scale); SciPy uses c (shape) and scale (\(=\beta\)).

• Numerical stability:

  • Use logpdf instead of pdf for likelihood work.

  • For extreme right tails, prefer sf/logsf instead of 1-cdf.

• MGF domain: the MGF exists only for sufficiently small positive \(t\); avoid relying on it outside its radius of convergence.

• Goodness-of-fit after fitting: KS p-values are not exact when parameters are estimated from the same sample; consider bootstrap.

12) Summary

• fatiguelife (Birnbaum–Saunders) is a continuous distribution on \((0,\infty)\) used for fatigue failure times.

• It has a clean normal connection via the transform \(\tfrac{1}{c}(\sqrt{x/\beta}-\sqrt{\beta/x})\).

• Mean/variance/skewness/kurtosis are available in closed form; the median equals \(\beta\).

• Sampling is easy with a NumPy-only transformation from a standard normal.

• In practice, scipy.stats.fatiguelife provides pdf, cdf, rvs, fit, and numerical entropy.